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Abstract. We present a one-equation subgrid scale model that evolves the turbulence energy corresponding to unresolved 
velocity fluctuations in large eddy simulations. The model is derived in the context of the Germano consistent decomposition of 
the hydrodynamical equations. The eddy- viscosity closure for the rate of energy transfer from resolved toward subgrid scales 
is localised by means of a dynamical procedure for the computation of the closure parameter. Therefore, the subgrid scale 
model applies to arbitrary flow geometry and evolution. For the treatment of microscopic viscous dissipation a semi-statistical 
approach is used, and the gradient-diffusion hypothesis is adopted for turbulent transport. A priori tests of the localised eddy- 
viscosity closure and the gradient-diffusion closure are made by analysing data from direct numerical simulations. As an a 
posteriori testing case, the large eddy simulation of thermonuclear combustion in forced isotropic turbulence is discussed. 
We intend the formulation of the subgrid scale model in this paper as a basis for more advanced applications in numerical 
simulations of complex astrophysical phenomena involving turbulence. 
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1. Introduction 

In the last decade, the significance of turbulence in various as- 
trophysical phenomena from stellar to cosmological scales has 
been recognised. In retrospect, this is hardly surprising, since 
virtually all the matter in the Universe is fluid, whereas the solid 
state is encountered as a rare exception. Moreover, both the in- 
tegral length L and the characteristic velocity V in astrophysical 
systems are quite large compared to terrestial standards, while 
the viscosity v is comparable to what is found for liquids or 
gas on the Earth. For this reason, the dimensionless Reynolds 
number 

Re = LVIv (1) 

becomes very large. A typical figure is Re ~ 10^"* for turbulent 
stellar convection. 

From the computational point of view, the number of de- 
grees of fre edom in a fluid dyriamica l system is given by the 
relation (see lLandau & LifshitJ 19871) 

~ (L/7?k)' ~ Re^/^ (2) 

The length scale ?7k is called the Kolmogorov scale and speci- 
fies the smallest dynamically relevant length scale. Due to the 
restriction of the CFL time step for compressible flow, the total 
number of operations required to compute the evolution over 
one sound crossing time is of the order 

A^ci- ~ A^"*^^ ~ Rel (3) 



Hence, N^^ ~ 10 operations would be required to solve the 
problem of stellar convection in a direct numerical simulation. 

The relevance of the Landau criterion |2l has been ques- 
tioned recently on the grounds of the intermittency of turbu- 
lence (Kritsuket al. 2005). In fact, the number of degrees of 
freedom refers to the ensemble of turbulent flow realisa- 
tions. At any particular instant of time, however, turbulent dy- 
namics is concentrated in regions of fractal dimension D < 3. 
Topologically, these regions can be either vortex filaments in 
subsonic flow or shocklets in the case of supersonic turbulence. 
The fractal dimension of vortex filaments is D = 1, whereas 
p - 2 for shocklets. According to the fi model (see iFriscl] 
1199.4 Sect. 8.5), the effective number of degrees of freedom is 
then given by 

A^eff ~Re'^^<^^". (4) 

If an adaptive numerical scheme with maximal efficiency in 
tracking the intermittent turbulent regions were applied, the to- 
tal computational cost could be lowered substantially in com- 
parison to a direct numerical simulation on a static grid. For 
example, it would appear feasible to treat subsonic turbulence 
up to a Reynolds number of 10"^ with an anelastic adaptive 
code on high-end platforms in the near future. On account of 
the limited efficiency of adaptive schemes, however, the ac- 
tual constraint might be lower. Apart from that, gravity, mag- 
netohydrodynamic effects and, possibly, reaction networks in- 
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crease the work load dramatically in astrophysical applications. 
Moreover, we have A^ci- ~ A^eifA^'^^ ~ Re"^"^ for the partic- 
ularly interesting case of supersonic turbulence. As a conse- 
quence, even with sophisticated adaptive schemes, it remains 
intractable to resolve completely the turbulent fluid dynamics 
encountered in astrophysics. 

In large eddy simulations (LES), on the other hand, only 
a limited number of degrees of freedom, which correspond 
to the largest scales of the system, is treated explicitly. For 
the turbulent dynamics on smaller scales, a so-called subgrid 
scale model is utilised. Among astrophysicists, the most of- 
ten used subgrid scale (SGS) model is numerical dissipation. 
This means that all fluctuations on length scales smaller than 
the resolution A of the numerical grid are smoothed out and 
it is assumed that the dynamics on length scales larger than A 
are more or less independent of the smaller scales. This point 
of view is m otivated by the second similarity hypothesis of 
iKolmogorovl U941) . which holds that the actual mechanism 
of dissipation is insignificant, provided that there is sufficient 
scale separation. In other words, on length scales Z ^ A, it is 
unimportant whether energy dissipation is caused by the mi- 
croscopic viscosity v at the length scale t/k or by numerical 
effects at the cutoff length A. The notion of numerical dissi- 
pation has been exhaustively investigat ed for the piece-wise 
para bolic method (PPM) proposed by IColella & WoodwardI 
f 1 984 ), which is one of the most popular finite-volume schemes 
applied in astrophysics. At least for statistically stationary 
isotropic turbulence, the numerical dissipation of the PPM ap- 
pears to work well as an implicit SGS model dSvtine et alJ 
boOO: SchmidL et al. 2005b). 

For the treatment of transient or inhomogeneous turbulent 
flow, however, an explicit SGS model becomes mandatory. One 
of the most prominent examples in contemporary theoretical 
astrophysics are numerical simulations of turbulent comb us- 
tion in type la supernovae ("Hillebra ndt & Niem ever"2000'). In 
this case, the velocity scale associated with SGS turbulence 
determines the eff'ective propa gation speed of flame fronts 
/Niemev er & HiUebrand3ll995l) . Part II of this paper will be 
dedicated to the problem of type la supernova simulations. 
The exchange of energy between resolved and subgrid scales 
is expected to become dynamically significant in the case of 
highly intermittent turbulence, for instance, in collapsirig tur- 
bulent gas clouds in the interstellar medium llLarsonll2003l 
iMac Low & K lessen'200?). 

In this paper (paper I), we present a general framework for 
the form ulation of SGS m odels based upon the, filtering ap- 
proach of'Germanol lll992h . The mathematical operation of fil- 
tering smoothes the flow on length scales smaller than the pre- 
scribed numerical resolution A. Consequently, a scale separa- 
tion is introduced, where the smoothed density, velocity, etc. 
are identified with the resolved quantities computed in LES. We 
have extended this formahsm to co mpressible flow s, using the 
Reynolds stress model proposed b v 'Canutol ll 1 997l) as a guide- 
line. N ext we discuss the one-equ ation SGS turbulence energy 
model llSchumannlll9"7i ISagauJEooii For the energy trans- 
fer across the numerical cutofl", we introduce a localised eddy- 
viscosity closu re which makes use of the dynamical procedures 
introduced bv iGermano et alJ lll99ll) . Hence, the SGS model 



becomes independent of a priori structural assumptions, in par- 
ticular, whether the resolved flow is homogeneous or stationary. 
However, it remains a necessary condition that turbulent re- 
gions become locally isotropic on length scales comparable to 
the numerical cutoff length A, because the linear eddy-viscosity 
closure presumes alignment between the turbulence stress and 
the rate-of-strain tensors. In the future, multi-parameter clo- 
sures for the turbulent energy transfer which are not subject 
to this restriction might be adapted. For the still more compli- 
cated non-local transport of kinetic energy on subgrid scales, 
we use a simple gradient-diffusion closure. In contrast to what 
has been suggested in the literature ( Pope ■2000.. Sect. 10.3), 
we find that the optimal diffusivity parameter is larger by about 
one order of magnitude than the SGS viscosity parameter, i.e. 
the turbulent kinetic Prandtl number is large compared to unity. 
Both the locaUsed eddy-viscosity and the statistical gradient- 
diffusion closure, respectively, were tested by means of data 
from simulations of forced compressible turbulence. As a case 
study, we present the LES of turbulent combustion in a periodic 
box. Although gravitational effects on subgrid scales, in prin- 
ciple, can be incorporated into the model as well, in paper I we 
restrict the detailed formulation and application to the case of 
negligible gravity. However, a simple closure which accounts 
for unresolved buoyancy effects in simulations of thermonu- 
clear supernovae will be discussed in paper II. 

2. Decomposition of the hydrodynamical 
equations 

Large eddy simulations pose the problem of scale separation. 
The numerically computed flow can conceptually be defined 
by a set of smoothed fields which correspond to the low-pass 
filtered physical flow realisation, where - n/A is the cutoff 
wavenumber for a numerical grid of resolution A. A low-pass 
filter is a convolution operator which is defined by 



q(x, t) - {q)G = I d jc G{x - X , t)q{x , f) 



(5) 



for a particular kernel G{x - x',t). The Fourier transform, the 
so-called transfer function G(k, f), drops to zero for wavenum- 
bers k ^ kc- Consequently, only modes of wavenumbers less 
than kc contribute significantly to the filtered field q{x, t). The 
exact, unfiltered variable q(x, t) corresponds to the limit k^^ 

oo 

oo. The filter operation smoothes out the fluctuations q — q-q 
on length scales I smaller than A. Albeit being mathematically 

OO 

determined by some dynamical equation, q{x, t) is generally 
not computable and therefore will be referred to as an ideal 
quantity. 

The dynamical equation for the ideal velocity field is the 
generalisation of the Nav ier-Stokes equation for compressible 
fluids (see .Landau & Lifshitz 1987i) : 



Deo - 

Df •' 



(6) 



The differential operator on the left hand side is the Lagrangian 
time derivate 



D 5 oo 
— = — + f ■ V, 
D? dt 



(7) 
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and the effective force density acting upon the fluid is defined 
by 

F^pf^-VP + V-o- + F'''^'\ (8) 



The first term on the right hand side is the pressure gradient, the 
second term accounts for viscous dissipation and the third term 
is the total external force per unit volume which encompasses 
gravitational and, possibly, stirring forces: 



The viscous dissipation tensor is defined by 

oo MOO°°^ 0O0o/'>=' 1°= 

o-ij = 2pvSij = 2pv 1 5 ij - -ddij 
where v is the microscopic viscosity of the fluid, 



Sij 



1 



oo „oo 

dui ovj 
dxj dxi 



(9) 



(10) 



(11) 



oo 

are the components of the rate-of-strain tensor and d - OiVi is 
the divergence of the ideal flow. 

Equation (|6} can be written in a conservative form as 



3 oooo oooo oo °° 

— pv + \-pv®v - F. 

at 



(12) 



The equality of the left hand side in both equations follows 
from the continuity equation which expresses the conservation 
of mass 

U CO oooo 

-p + V-pv^Q (13) 
of 

The goal of the filtering approach is the formulation of dy- 
namical equations for smoothed quantities which are amenable 
to the numerical computation. We define the Favre or mass- 
weighted filtered velocity field 



v{x, t) = 



OOOO 

{pv}c 



{p)c (14) 

= r d\'G{\x-x'\)^{x\t)v{x',t), 

p(x, t) J 

CO 

where p — {p)g is the filtered mass density. For a homoge- 
neous and time-independent kernel, the filter operation com- 
mutes with the Lagrangian derivative. Favre filtering the dy- 
namical equation Mil , one obtains 



3 oooo oo 

—pv + V-{pv®v)G ^F, 
3t 



(15) 



where F - {F)g- The filtered equation can be cast into a form 
analogous to equation (|6|l, 

OOCO CO 

p—v^F + '^-T{pv,v), (16) 
by virtue of the generalised turbulence stress tensor 

oooo oo oooo oo 

T{pv,v) - -{pv ® v)c + pv®v (17) 



which lOermanol ( 1 1 992l) introduced without mass-weighing for 

COCO CO 

incompressible turbulence. We will use = T(py,-, vj) as a 
shorthand notation for the components of the turbulence stress 
tensor. The term 3jTij in the momentum equation stems from 
the non-linear advection term in the Lagrangian time derivate 
and can be interpreted as the stress exerted by turbulent velocity 
fluctuations smoothed out by the filter. In the following, we 
will encounter a variety of r-terms. For this reason, we define 
t(-, ■) as a generic bilinear form which maps any pair of ideal 
fields to a mass-weighted smoothed field which is called the 
generalised second moment. The resulting field can be scalar, 
vectorial or tensorial. Of course, the notion of a generalised 
moment applies to moments of higher order as well. 

oo 

A dynamical equation for the specific kinetic energy, k - 
^\v\ , is readily obtained from the contraction of equation (|6) 

oo 

with the ideal velocity field v: 
D 



p — k = V F. 
^Dt 



(18) 



The mass-weighted filtered kinetic energy k(x, t) is defined by 



k{x, t) - 



{pk}c 

CO 

{p)c 



(19) 



Filtering equation J18> results in the following equation for 

k{x, t): 



p—k - V ■ ;f*'"^ = d ■ f + V ■ 

Df 



OOCO CO 



V ■ rip v,v) 



(20) 



In addition to the turbulence stress term on the right hand side, 
there is a non-local transport term which is given by the diver- 
gence of the turbulent kinetic flux 'J^^^^'^\ The flux is defined by 
the contraction of the completely symmetric generalised third- 

oooo oo oo 

order moment r(p«, «, «). In component notation, we have 



2Tl 



(kin) 



OOOO oo oo 



Tijj = TipVi, Vj, Vj) 

OOOO OO OO oooo CO 

-{pViVjVj}c + {pViVj)cVj ■ 



(21) 



Since the filtered kinetic energy is a second-order moment 
of the ideal velocity field, contributions from velocity fluctua- 
tions on all scales are included. For this reason, k differs from 
the specific kinetic energy of the smoothed flow, jl^P, and 



(22) 



can be identified with the generalised turbulence energy asso- 
ciated with scales smaller than the characteristic length of the 
filter ( )q. The dynamical equation for A^,urb is obtained by sub- 
tracting 



D /I 



p — -\v\- \-v 



F + V ■ T(pv, v) 



(23) 



from equation OQI . The result, in component notation, reads 



,ooco oo 



P^^turb - 5,^F.*'"> = Tipvi, vj)S ij -r{ vi, Fi), (24) 
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where t(i',, Fj) is the contraction of the tensor 

00°° 00 °° 
t(v,F) - -{v ® F)c + v®F. 



and the transport term D is given by 



(25) 



In order to put equation (I24> into a form which is more ad- 
equate for modelHng purposes, flux terms are spHt off on the 
right hand side. Let us first substitute the definition of the ef- 
fective force 



cx) „ 00 



t{v,, Fi) = -T(Vi, diP) + T(Vi, dj(T,j) + T(Vi,F''^'>). (26) 

The three resulting generahsed moments respectively corre- 
spond to pressure, viscous and external forces. Because stir- 
ring forces usually supply energy on the integral length L of 
the flow only, it follows that T(t;,, F.'^*) is negligibly small for 
A <K L. The first and the second term on the right hand side of 
equation (I26> can be split as follows: 



00 CX) 00 



-T(Vi,diP) = -diT(Vi,P) + T(d,P), (27) 

00 00 CO CO CO 

T(vi,djcrij) = djT( vi,crij) - T{Sij,aij). (28) 



The new r-terms on the right hand side are, respectively, 

the pressure-dilatation and the rate of viscous dissipation. 

Substituting the expression for the viscous dissipation ten- 
sor ( I10> . it follows that 



CO coco coco 

TiSij,o-ij) = -{pv\S I }G-2{pvSij}cS,j. 



(29) 



The norm \S*\ is defined by the total contraction of the trace- 
free part of the rate of strain tensor: 



CO CO 00 00 

\s'\'^2s;js;j^\s\'--d\ 



(30) 



The norm of the rate-of-strain tensor is a probe of the velocity 
fluctuations on the smallest dynamical length scales which are 
of the order of the Kolmogorov scale t/k- We shall assume that 
the characteristic length of the filter is much greater than the 
Kolmogorov scale. In this case, the rate of strain of the filtered 
velocity field is much less than the rate of strain of the ideal 

velocity field, i.e. \S*\^ ^ Consequently, the first term on 
the right hand side of equation ( I29t dominates the second term, 
and we can set 

T(Sij,crij)^-{yp\S'\^)c. (31) 
Summarising, equation J24> can be written in the form 



D 

P^^tuib - 2) = r + I - -H e), 
where the source contributions on the right hand side are 



,00 CO CO , 



r = <p, Vigi)c - pvigi, 

2 - TijS ij, 

COCO 

pA = -{dP)c + dP, 

MOO ^„ 1 

pe^{vp\S*\^)c, 



(32) 



(33) 
(34) 

(35) 

(36) 



dxi 



1 00 00 



The generalised moment 

Al = -{vP)g + vP 



(37) 



(38) 



accounts for the transport of turbulence energy due to pressure 
fluctuations. 

For the internal energy, the filtered dynamical equation is 



D 



cooo 00 



(39) 



^PXVT + Q-{Pd)c-pe. 



The source term Q accounts for heat generation by chemical or 
nuclear reactions. The transport of heat due to turbulent tem- 
perature fluctuations gives rise to the generalised conductive 
flux, 

l^^„A\ COCO OOCo'^ 

CO 

for fluid of thermal conductivity x- The additional transport 
term on the left hand of equation J39t side arises from the trans- 
port of heat by turbulent advection. In the case of buoyancy- 
driven turbulence, this transport mechanism is known as con- 
vection. The generalised convective flux is defined by 



r'""'' ^Tipv,h) = Tipv,ei„,)+fi, 



(41) 



where h - + P/pis the specific enthalpy. 

Adding the budgets of the specific kinetic energy and 
internal energy eim, we obtain the total energy per unit mass 
on length scales I ^ A, i.e. etot = ^int + \\v\'^- The dynamical 
equation governing the evolution of etot is 

D 



V ■ |^y^(':o"'i) <p^(conv) _ ^pj 

= p;rvr + e+p(^ + e) 

/ 0000 00 

g + f'> + V-T(pV,v) 



(42) 



+ V ■ 



This conservation law in combination with equation i32i ex- 
tends the Gennano consistent decomposition to compressible 
turbulence. The sum of internal energy and kinetic energy on 
all scales is etot + ^tmb = ^int + k. Comparing equations ( I32t and 
(I42> . one can see that p{A + e) accounts for the dissipation of 
turbulence energy into internal energy by compression effects 
and viscous dissipation, respectively. The turbulence produc- 
tion term E is related to the energy transfer through the turbu- 
lence cascade across the characteristic length of the filter The 
injection of energy due to buoyancy and the action of stirring 
forces on length scales larger than A is given hy v ■ {g + /*^*], 
whereas the interaction of gravitational potential energy fluc- 
tuations and turbulence on length scales smaller than A is ac- 
counted for by the term F. 

From the discussion in this Section, it should become clear 
that the presumed scale separation in LES is essentially based 
upon the disentanglement of a variety of dynamical effects. 
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This task is considerably complicated by the non-linear inter- 
actions across the cutoff scale, which become manifest in the 
various generalised moments occurring in the equations i20\ . 
i3H . i39\ and J42l i. Hence, one faces the problem of finding 
closures for the generalised moments in the decomposed dy- 
namical equations. 

In the simplest of all cases, a closure is a sufficiently con- 
vincing argument for neglecting a certain term. This kind of 
closure is applied in many cases. In the proper sense, a clo- 
sure is a more or less tentative approximation which is made on 
grounds of heuristic physical arguments. Two major categories 
of closures can be distinguished: An algebraic closure is some 
function of filtered quantities only. Usually, algebraic closures 
contain at least one free parameter Depending on whether this 
parameter is a constant or varies in space and time, the clo- 
sure is either statistical or locaUsed. On the other hand, dy- 
namical closures determine generalised moments from addi- 
tional dynamical equations. However, these equations, in turn, 
entail closures for higher-order generalised moments. This is 
the problem of the infinite hierarchy of equations for filtered 
quantities. Inevitably, the hierarchy must be truncated at some 
point with the help of algebraic closures. 

3. The subgrid scale turbulence energy model 

In LES, numerical solutions for the filtered quantities p, P, T, 
D and e,o, have to be computed from the continuity equation 
for p, the momentum equation il6i . the energy conservation 
law i42i and the equation of state. The filter operation naturally 
introduces a cutoff which is related to the numerical resolution. 
Since the filtering in LES is not necessarily explicit but some- 
times inherent to the numerical scheme, we will subsequently 
use the generic notation ( )eif . For example. 



v(x, t) = 



{pv{x,t))^s 

CO ' 

<P>eff 



(43) 



is the mass-weighted velocity field which is to be computed nu- 
merically. For finite precision, the numerical solution actually 
corresponds to a whole ensemble of exact flow realisations. In 
this regard, one can think of a reduction of the number of de- 
grees of freedom due to filtering. 

The length scales smaller than the characteristic length Agfi 
of the effective filter are the subgrid scales (SGS). The scales 
I <: Aeff, on the other hand, are numerically resolved. A com- 
plete set of closures for the generalised moments in the dy- 
namical equations constitutes the subgrid scale model. A gen- 
eral SGS model which includes dynamical equations for the 
moments of second and third order was formulated by Canuto.. 

Unfortunately, the computational cost of solving the 
whole set of dynamical equations for the generalised moments 
is considerable. Moreover, the problem of stability appears to 
be non-trivial. 

The SGS model which will be discussed in the following 
involves the solution of the dynamical equation for the subgrid 
scale turbulence energy only. The definition of the density of 
SGS turbulence energy is as follows: 



_ 1 2 



1 



The magnitude of SGS velocity fluctuations is given by q^g^. 
The equation governing the evolution of the specific turbulence 
energy ^sgs = K^^^jp is just equation j32t with the filter ( )eff. 
For the various second-order moments in this equation we wiU 
invoke standard algebraic closures. Hence, the implementation 
of the SGS model requires the solution of only one additional 
dynamical equation. The inherent hmitations of the simple clo- 
sures are in part compensated by localisation. Thereby, the SGS 
model becomes basically independent of a priori model pa- 
rameters which presume certain flow properties such as homo- 
geneity. In essence, the SGS model which will be formulated 
is robust and particularly well suited for complex flow geome- 
tries and transients, although requiring relatively little compu- 
tational resources. 

There are two different sources of SGS turbulence produc- 
tion. The first one is the SGS gravity term Fsgs ( I33> which ac- 
counts for the conversion of potential into kinetic energy and 
vice versa due to correlations between SGS fluctuations of the 
velocity and gravity. A putative closure for SGS buoyancy in 
reactive flows will be presented in paper II. The other produc- 
tion term is the rate of energy transfer Esgs - TijS ij across 
the length scale Agfi due to non-linear turbulent interactions. 
In general, energy transfer is the primary source of SGS tur- 
bulence. A common closure is based on the eddy-viscosity hy- 
pothesis f or the the tra ce-free part of the SGS turbulence stress 
tensor (cf. lPopel2000l Sect. 10.1.): 



where 



T-j = 2pvsgs5^ = 2pv,o,\Sij - ^d6ij 



1 2 



(45) 



(46) 



3 ' ' 3 
This closure is formulated analogously to the viscous stress 
tensor in a Newtonian fluid. The eddy viscosity Vsgs is assumed 
to be proportional to the product of the characteristic length 
Aetf of the numerica l scheme and the characteristi c velocity of 
SGS turbulence (cf. ISagautlEoOll Sect. 4.3, and IPonelEoOol 
Sect. 13.6.3), i. e.. 



(47) 



sgs — -i"isgs 

The length scale (y - CyAeff / V2 is thus associated with SGS 
turbulence production. 

Among the dissipation terms, the rate of viscous dissipa- 
tion esgs defined in equation ( I36> dominates in subsonic turbu- 
lent flows. Assuming a Kolmogorov spectrum, the mean SGS 
turbulence energy corresponding to a sharp spectral cut-off can 
be related to the mean rate of dissipation: 

^^£(A:)dA:= -C<esgs)'^'(^) ■ (48) 
Hence, assuming fliat C ~ 1.65 jYeung & Zhoul 19971) . 



= 2 — ^ -0.81 ^ . (49) 

\ 2 / A A 



(^sgs) 



Conjecturing that the above relation also holds locally (cf.lPoi 
.2000i Sect. 13.6.3), we have 



{p\W%«-p\v[ 



(44) 



''sgs H: 



sgs 



Aeff 4 



(50) 
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where — 2 V2AefF/C(; and Q ~ 1. Basically, equation (|50| 



implies that SGS eddies of kinetic energy ~ q^^^ are dissipated 
on a time scale ~ ^e/^sgs- 

Pressure dilatation poses severe difficulties because one 
needs to find the correlations between pressure fluctuations and 
compression or rarefaction of the fluid. The first-order hypoth- 
esis is that kinetic energy is dissipated if the fluid is contracting 
(d < 0). In the opposite case (d > 0), internal energy is con- 
verted into mechanical energy which produces turbul ence. This 
line of reasoning leads to the closure proposed by IPeardor^ 
lll973h : 

/Isgs = C^^sgs^^- (51) 

Unfortunately, numerical tests reveal that this closure is ex- 
tremely crude. Since compressibility eff'ects tend to diminish 
toward smaller length scales, the above closure will do for the 
LES of subsonic turbulence. In the case of supersonic turbu- 
lence, however, A^gs becomes more significar it. Alternative c lo- 
sures for pressure-dilatation are described in 'Canutol (ll997h . 

A customary algebraic closure for the transport term in 
equat ion ( I32t is the gradient-diffusion hypothesis (cf. ISagautI 
1200 1[ Sect. 4.3)1 



d_ 



dk^g 
dxi 



2 dq^g^ 



(52) 



The characteristic length scale of diffusion is defined by t,^ = 
C/fAeff/ V2 and the SGS diffusivity is given by Ksgs - ^«g'sgs. 
The notion of a turbulent diffusivity of kinetic energy stems 
from the analogy to the thermal diffusion of internal energy on 
microscopic scales. This analogy also suggests the definition of 
a kinetic Prandtl number. 



Vsgs Cv 

Pridn = = — ■ 



(53) 



Summarising, if gravitational effects on subgrid scales are 
negligible, we obtain the following dynamical equation for the 
SGS turbulence energy: 



(54) 



The remaining problem is the determination of the closure pa- 
rameters Ck, Cv, Ca and which are dimensionless similarity 
parameters, i.e. the values become asymptotically scale invari- 
ant in statistically stationary isotropic turbulence. 

4. Closure parameters 

In this section, methods for the calculation of closure parame- 
ters will be discussed. In particular, we will present a so-called 
dynamical procedure for the computation of the eddy-viscosity 
parameter Cv Originally introduced by engineers in order to 
improve the performance of simple SGS models such as the 
Smagorinsky model for turbulent flows near walls, the appli- 
cation of dynamical procedures for the localised computation 

' Also known as Kolmogorov-Prandtl relation. 



of closure parameters has turned out to be a powerful tool for 
the treatment of inhomogeneous and non-stationary turbule nce. 
For this reason, we adapted a procedure proposed bv Kim et alJ 
( 199% for the LES of turbulent combustor flows. Using data 
from simulations of forced isotropic turbulence, we found that 
this procedure yields a significantly better match with the rate 
of production than the statistical closure with a constant param- 
eter. For the parameter of dissipation, Q, we propose a semi- 
statistical solution: A time-dependent value is determined from 
the energy budget of the resolved flow in extended spatial re- 
gions. Regarding the non-local transport, the gradient-diffusion 
closure produces satisfactory results if the parameter C^ is de- 
termined appropriately. 



4.1. Production 

The rate of production Esgs corresponds to dissipation of ki- 
netic energy on resolved scales due to the effect of subgrid 
scale turbulence. Pictorially, unresolved eddies drain energy 
from larger eddies at the rate Esgs. This idea motivated the eddy- 
viscosity closure (I45> for Ssgs- Extending further the analogy 
between viscous and turbulent dissipation, an experimental as- 
sertion known as the law of finite dissipation could be carried 
over to the production of SGS turbulence: If, in a large eddy 
simulation of turbulent flow, all the control parameters are kept 
the same except for Aeff , which is lowered as much as possible, 
the energy dissipation per unit mass, Ssgs, behaves in a way 
consistent with a finite positive limit^. This suggests that the 
parameter Cy in the definition of the turbulent viscosity ( I47> 
becomes asymptotically scale-invariant in the limit Agfj /L — > 
and assumes a universal value in the stationary limit. 

We verified this hypothesis by analysing data from nu- 
merical simulations of forced isotropic turbulence. The driv- 
ing force which supplies energy on the characteristic length 
scale L is modelled by a stochastic pro cess in Fourier space 
dEswaran & PoDelll988l: rSchmidlll2004 . Under the action of 
this force, the flow evolves on the characteristic time T which 
is called the large-eddy time scale. In the statistically stationary 
limit, the flow velocity is of the order V - LjT .\n addition, the 
weight of solenoidal (divergence-free) relative to dilatational 
(rotation-free) components of the force field can be varied by 
setting the control parameter ^ in the range between 1 and 0. 
Choosing different characteristic Mach numbers V/cq, where cq 
is the initial sound speed, and values of ^, we performed sev- 
eral simulations using the piece-wise parabolic metho d (PPM) 
with N = 432^ grid cells dColella & Woodward 19841) . A real- 
istic equation of state f or electron-dege nerate matter was used 
in these simulations (see lReuieckel200ll) and the numerical dis- 
sipation of PPM provided an implicit subgrid scale model. 

It is possible to evaluate generalised moments from the sim- 
ulation data on a length scale which is large compared to the 
cutoff scale A. To that end, let us introduce new smoothed fields 

and which are associated with a basis filter ( )< of char- 



- The formulation is the same as in lFrischlil99l . beginning of Sect. 
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Table 1. Closure parameters for selected flow realisations from 
three different simulations of forced isotropic turbulence 



V/co 




tIT 


A</A 


<c,,> 






Prkin 


0.42 


1.0 


2.5 


6.9 


0.0512 


0.0653 


0.401 


7.8 


0.66 


0.75 


2.0 


9.6 


0.0447 


0.0611 


0.401 


9.0 


0.66 


0.75 


4.0 


6.8 


0.0476 


0.0649 


0.422 


8.7 


0.66 


0.75 


9.0 


6.5 


0.0451 


0.0661 


0.529 


11.7 


1.39 


0.2 


3.5 


16.1 


0.0370 


0.0986 


0.481 


13.0 


1.39 


0.2 


6.0 


6.4 


0.0422 


0.0806 


0.515 


12.2 


acteristic lenj 


;thA<: 


<P><, 


p<v< 


oooo 

= {pv)<. 




(55) 



If A< is sufficiently large compared to Agff, then q = {q}< ^ 
{q)t: for any quantity q. This property of low-pass filters be- 
comes immediately apparent in spectral space in which the fil- 
ter operation is just a multiplication of Fourier modes with the 
transfer function. Consequently, the turbulence stress tensor at 
the level of the basis filter is approximately given by 



T<(pi;, d) k. T^ipv, v) - -{pv ® d)< + p V ® V 



(56) 



Evaluating T<(pv,v), it is possible to determine the eddy- 
viscosity parameter from the rate of energy transfer across the 
length scale A< : 



c: 



(57) 



where 2< = T%(pui,Vj)Sfj and p^k^^^^^ ^ -ir<(p!;;,i;,). We 
computed Cy from several flow realisations. A sample of av- 
erage values for the whole domain is listed in Table ^ We 
found {Cy} * 0.05 for dev eloped turbulence, in agreement with 
the hterature llPopell200(il) . Only in the case of transonic flow 
with a mostly compressive driving force the eddy-viscosity pa- 
rameters appears to be systematically lower Fig.n(a) shows 
a visualisation of the turbulence energy isosurfaces given by 
^turb ~ 0.25V^ with the contour sections of the dimensionless 
rate of energy transfer £< = (T/poV^)!.^ for V/cq = 0.66 at 
time t = 4T. The corresponding contours obtained with the 
eddy-viscosity closure and Cy - 0.0476 are plotted in panel 
(b). Clearly, the rate of energy transfer is not well reproduced. 
Although there is a significant correlation of about 0.8, the 
magnitude of spatial variations is greatly reduced. 

In fact, Cy exhibits spatiotemporal variations comparable 
to the mean value. In consequence, the assumption of a constant 
eddy-viscosity closure parameter is not valid. However, the in- 
formation about the variation of Cy is not available in a LES. 
A solution to this problem can be found by means of a similar- 
ity hypothesis which relates the energy transfer across different 
scales. Let us consider a length scale Aj which is somewhat 
larger than A<. Introducing a suitably defined filter operation 
( )t of characteristic length Aj, the turbulence stress Tjipv, v) 
is given by an expression analogous to the right hand side of 
equation (l56l : 



Here p^^' - (p"^)t and d*^* = (p*i'*)t/(p^)t- The stress tensors 
associated with the length scales Aj an d A<, respective ly, are 
related by an identity which Germano (iGermanoll 19921) origi- 
nally formulated for incompressible turbulence: 



Tj(pv,v) = {tJpv,v)}t: + tt:(p^v^,v''). 



(59) 



The first term on the right-hand side is the filtered turbulence 
stress tensor associated with the length scale A<, whereas the 
second term accounts for the turbulence stress on intermedi- 
ate length scales in between A< and Aj. For small scaling ra- 
tios yj - At/A<, there is significant correlation not only be- 
tween Tjipv, v) and (r^ipv, v)}j, but also between Tj ipv, v) and 
TTip^v "^, v^)- In particular, this was demonstrated bv lLiu et"al] 

(Il994{) from velocity measurements in round jets . 

I i . : I upon these experimental findings. .Kim et alJ (Il999l) 
proposed a similarity hypothesis for the eddy-viscosity param- 
eter: 

y(T) 

(60) 



pmATkl'^\S*m\' 



where |5 **^^| is the norm of 



(61) 



,(T)„(T) ^ j,(T) 



(58) 



and S*^^' = T^(fi^uf, v'j)S^^\ The specific turbulence energy kj 
corresponding to intermediate velocity fluctuations in between 
the basis and the test filter level is defined by 

p(T)^^ = -^TT(p<vf, Vf) = -^TTipVi, Vi) - <p"Cb>T- (62) 

The second expression for kj follows from the contraction of 
the Germano identity ( I59> . Thus, the parameter Cy for the 
eddy-viscosity closure at the level of the basis filter is deter- 
mined by probing the flow at the length scale Aj > A<. This is 
why ( )t is called a test filter. 

Using data from the simulations of forced isotropic turbu- 
lence, we tested the proposition made above by computing ex- 
plicitly the rate of energy transfer across a certain length scale 
A< and comparing it to the eddy-viscosity closure with the clo- 
sure parameter calculated at test filter levels for different scal- 
ing ratios yx- In order to apply approximation i56\ . we had 
to choose a basis filter length A< which was at least an or- 
der of magnitude larger than the resolution A in the simula- 
tions. On the other hand, a sufficient range of inertial length 
scales greater than A< is required for the test filter operation. 
These requirements substantially constrained the choice of A<. 
Further complications come from the so-called bottleneck ef- 
fect which causes a distortion of the energy spectrum function 
for wave numbers close to the cutoff at k^ - n/A (Dobler et al-j 
l2003UHaugen & BrandenburJ2004 . A detailed discussion of 
the kinetic energy spectrum functions and, particularly, the bot- 
tleneck effect in turbu lence simulations with PPM is given in 
ISchmidt etaD(l2005bl) . As one can see in Fig.|2 the match be- 
tween the probability density functions of the dimensionless 
rate of energy transfer and the corresponding localised closure 
is substantially better than for the closure with constant eddy- 
viscosity parameter This is highlighted by the statistical mo- 
ments listed in Table |2l In particular, the variance of the en- 
ergy transfer is largely underestimated by the statistical closure. 
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fig.transf .png fig_closr_ave.png 
(a) Explicit (b) Statistical closure 



f ig_closr_locl . png £ig_closr_locl_nb . png 

(c) Localised closure (d) Localised closure excluding backscattering 

Fig. 1. Isosurfaces of the turbulence energy k^^^^^ - Q.25V^ with contour sections of the dimensionless rate of energy transfer. The 
flow reahsation is taken from a simulation with characteristic Mach number V/cq = 0.66 at time t = 4T. 



This also becomes apparent from the three-dimensional visual- 
isations in Fig.n(b) and (c), respectively, which suggest that 
variations of the energy transfer are flattened by a wide mar- 
gin in the case of a constant eddy-viscosity parameter, while 
the localised closure reproduces local extrema quite well. On 
the other hand, it appears that the characteristic length At of 



the test filter should not be chosen too large in relation to . 
Otherwise the mean of the energy transfer will be systemati- 
cally underestimated (Fig.|2and TableO. 

The variability of Cy is illustrated by the probability den- 
sity functions plotted in Fig.|3 The similarity of the functions 
suggest a fairly robust behaviour of Cj,^' for driven isotropic 



W. Schmidt et al.: A localised subgrid scale model for fluid dynamical simulations I 



9 




Dimensionless rate of energy transfer DImensIonless rate of energy transfer 



Fig. 2. Probability density functions for the rate of energy transfer £< across the length scale A< in dimensionless scaling at two 
different instants of time. 



Table 2. Statistical moments of the dimensionless rate of en- 
ergy transfer for the probability density functions plotted in 
Fig.El(b). 
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= 4.0) 


0.183 
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turbulence. In a fraction of roughly 15 to 20 % of the domain, 
negative values of the closure parameter are found which are 
commonly interpreted as inverse energy transfer from length 
scales smaller than Aj toward larger scales. This phenomenon, 
which is also know as backscattering, is predicted by turbu- 
lence theory. However, as we shall argue in Sect.|5] backscatter- 
ing introduces numerical difficulties in combination with PPM. 
But panel (d) in Fig. demonstrates that the localised closure 
is superior even when negative values of the eddy-viscosity are 
suppressed. 

For the application in LES, the basis filter corresponds to 
the effective filter introduced in the previous Sect., and the test 
filter is applied to the computed fields p{x, t) and v{x, t). Then 
we have 



(63) 



The characteristic length scale of SGS turbulence production, 
{y depends on the scaling ratio yj - Ax/Aeft and, con- 
sequently, dyCyAen / V2 is proportional to the parameter /3 - 
A^ff/A. ISchmidt et al. (2005b) determined /3 ~ 1.6 for the sta- 
tistically stationary turbulent regime in simulations with PPM. 




0.20 



Fig. 3. Probability density functions for the localised eddy vis- 
cosity parameter calculated from different flow reaUsations. 



4.2. Dissipation 

The localised closure for the rate of production works be- 
cause the energy transfer across a certain cutoff wavenumber 
is mostly determined by interactions between Fourier modes 
within a narrow band around the cutoff. Concerning the rate of 
dissipation fsgs, we encounter an entirely different problem. In 
fact, viscous dissipation takes place on length scales which are 
of the order of the Kolmogorov scale t/k ^ Ae^ . There is no ob- 
vious similarity between the dissipation on resolved scales (due 
to SGS turbulence) and the dissipation on subgrid scales (due to 
microscopic viscosity). The simplest of all SGS models, which 
is known as the Smagorinsky model, assumes a local equilib- 
rium between the dissipation on resolved and subgrid scales, 
respectively. However, it is the very point of the SGS turbu- 
lence energy model that such a balance does not hold locally. 
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Nevertheless, the mean rate of energy transfer can be related to 
the rate of viscous dissipation in the case of homogeneous tur- 
bulence. If the flow is inhomogeneous, equilibrium might be 
assumed at least for some nearly homogeneous regions. Thus, 
we attempt to determine the closure parameter from the av- 
eraged energy budget on the test filter level for a suitably cho- 
sen flow region. 

The m ethod is loosely based on the variational approach of 
ICihosal et al. ( 1995) . They subtracted the test-filtered SGS tur- 
bulence energy equation (I54t from the corresponding equation 
for the turbulence energy at the level of the test filter in order 
to determine as a function of both space and time. Our ap- 
proach is an intermediate one, where spatially averages energy 
equations are considered. For the mean SGS turbulence energy, 
averaging equation (|53 yields 



Due to the large number of filtered quantities, the complete nu- 
merical computation of the source terms in the above equation 
would be rather demanding. For this reason, we drop the con- 
tributions (11) and (111) while only retaining (I), which is pre- 
sumably the most significant contribution to the rate of energy 
transfer across Ax. Then {Kj} is approximately given by 



df 



(70) 



Invoking the closure dimensional closure i50\ both for the 
SGS rate of dissipation fsgs and the rate of dissipation at the test 
filter level, we obtain the following expression for the mean rate 
of dissipation on length scales in between Aetf and At: 



Here it is assumed that the surface contributions from the trans 
port term cancel out or at least can be neglected. Since 

d_ . \ I d 



D 



we also neglect the effect of advection by the resolved flow. 
The turbulence energy density associated with the characteris- 
tic scale of the test filter is defined by the trace of the Germano 
identity (|59} 

1 COOO CO 1 1 

--TT(pvi, vi) = --<t/,)t + -Tripvi, Vi) = {K^g^h + Kt, (66) 

where Kj - p'-^^kj i62i . The spatial average of the turbulence 
energy (I66> is given by the following dynamical equation: 



O / coco CO ITW 

-{K,g,+KT}^{Tj(pVi,V,)S^^^) 



- (p('tsgs + fsgs) +P*'^'(^T + ej)) ■ 



(67) 



Equations ( I64> and (I66> in combination with the Germano 
identity (|59j imply 



- {p^^HAt + ex)) . 



(68) 



Substituting the eddy-viscosity closures for the various produc- 
tion terms on the right-hand side, the above equation becomes 



(I) 

- <p<T)^x> + <P<^'eT> 

+ (<pVsgs^pT^,; ' ^^-pV.sg.sl^f} 
(U) 

- ^ (</:sgs)T«?<^> - K,g.d) . 

(Ill) 



(69) 



Furthermore, setting 



(72) 



the closure parameter is determined by 

[at ' \ ■ 

,3/2 



df 

xAt/p<t> 



A</rT> - (c,pm AT^/^l^mp) + (1 + c,| (^t^^^ 

<P^sgs>T 

-h/CTi -TTP^sgs) • 

(73) 



l(T) 



Contrary to the eddy-viscosity parameter Cy which varies both 
in space and time is a time-dependent constant for a suit- 
ably chosen spatial region. For homogeneous turbulence, there 
is only one region encompassing the whole domain of the flow. 
In a stratified medium, it is appropriate to average horizontally. 
Then varies with depth. For turbulent combustion problems, 
such as type la supernova explosions, one can distinguish fuel, 
the burning zone and the burned material within. For each of 
these three regions a value of the dissipation parameter is calcu- 
lated as a function of time. For the pressure-dilatation parame- 
ter Ca, on the other hand, we preliminarily assume the constant, 
time-independent v alue Ca - for subsonic turbulence (see 
iFurebv et alJl997l) . 

4.3. Diffusion 

As in the case of the energy transfer, we shall first consider the 
problem of non-local transport at the level of a basis filter of 
characteristic length A< which is large compared to the numer- 
ical cutoff length. The generalised kinetic flux (12 1> is given by 



r: 



(kin) < 



+ {pUiVj}^v^ -p^vfu'ju^ 



and the pressure diffusion fluxl38lreads 



(74) 



(75) 
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Assuming that the total flux vector 
the turbulence energy gradient Vfcj^j,^ 
closure can be written as follows: 



+ //"^ is aligned with 
the gradient-diffusion 



'"turb ^"turb- 



(76) 



Contracting the above relation with Vfej^^.^^ and averaging over 
the domain of the flow, one obtains 



/-'(vec)< _ 2_ 'u™ 



(77) 



A sample of values for C)^ is listed in Table[2 In agree- 
ment with a turbulent kinetic Prandtl number of the order unity, 
j^(vec)< the same ord er of magnitude as the closure param- 
eter Cy (see Pope 2000, Sect. 10.3). Contour sections of the 
flux magnitude |^(*"'^'^ + fi^\ and the coiTesponding closure at 
the k^^^^^ = 0.25 isosurfaces for V/cq = 0.66 at time f = 47" are 
shown in Fig. |3 However, as one can see from a comparison 
of the panels (a) and (b), the closure underestimates the diffu- 
sive flux by about an order of magnitude. Even more clearly, 
this is demonstrated by the probability distribution functions 
plotted in Fig. |5] We also investigated the hypothesis of set- 
ting the turbulent diffusivity p arameter equal to the localised 
eddy-viscosity parameter (see Sagaut 2001: Kim e t al.l [19991 
Sect. 4.3). Since negative diffusivity would induce numerical 
instability, we truncated the diffusivity parameter at zero, i.e. 

- Ct^'^. The resulting visualisation in panel (c) of Fig.|3 
and the corresponding graph in Fig.|5l however, show very little 
if any improvement compared to the statistical closure. 

The reason for the discrepancies is the flawed assumption 
of alignment between the turbulent flux vector and the energy 
gradient. Setting 



c 



(scl)< 



Vk 



* I Ik^ 
turb' -y '^^turb 



(78) 



where an equality of the flux magnitude but not the direction 
is presumed, results in significantly larger turbulent diffusivity 
(see Table [0. In particular, panel (d) in Fig. |4] and the proba- 
bility distribution functions shown in Fig.|5]reveal a very close 
match between the explicitly evaluated turbulent flux and the 
gradient-diffusion closure with the parameter C['''^'** = 0.422. 
Remarkably, the implied turbulent kinetic Prandtl number is of 
the order of ten rather than unity. 

It appears that the gradient-diffusion closure provides a dif- 
fusive mechanism which accounts for the intensity of turbulent 
transport but fails to reproduce anisotropic properties of third- 
order generalised moment. This is why advanced statistical the- 
ories of turbulence abandon the gradient-diffusion closure and 
introduce dynamical equations for the third-ord er moments or 
make use of other, more sophisticated closures dCanutol l 19971 
[(Tanuto & Dubovikov 1998). Such equations have be en sug- 
gested for the application in SGS models as well dCanutol 
11994 . On account of the difficulties solving these equations, 
however, we prefer the simple algebraic closure i52i with a 
constant diffusivity parameter 



0.8 



0.0 









- / / 

- / / 

" / / 
1/ / 




^^^^^^ 

- 

V/co = 0.66, <■ = 0.75 






t/T = 2.0, Ag/A = 1 6 


l' ^ 




explicit 






scolor closure - 






vectoriol closure 


i/ y 




(Tk = 1 (7t = 2-0) - 










0.00 



0.02 0.04 0.06 

DImensionless flux magnitude 



0.08 




V/co = 0.66, f = 0.75 
l/T = 4.0, Ag/A = 1 1 



explicit 

scolor closure 
vectoriol closure 
ay = 1 (7T = 2.0) 



0.00 



).02 0.04 0.06 

DImensionless flux magnitude 



0.08 



Fig. 5. Probability distribution functions for + and 

the coiTesponding gradient-diffusion closures with different 
turbulent diffusivity parameters. 



corresponding to the turbulent diffusivity 



(80) 



C. ^ 0.4 



(79) 



According to our numerical investigation, Q ~ 0.4 is represen- 
tative for stationary isotropic turbulence of Mach number ^ 1 . 
In the case of developing turbulence, the effects of turbulent 
transport are rather marginal, and the deviations introduced by 
the statistical diffusivity J80> are not overly important for the 
subgrid scale dynamics. For higher Mach numbers, however, 
there appears to be a trend towards systematically larger diffu- 
sivity. 

In a similar fashion as the gradient-diffusion hypothesis, a 
turbulent conductivity ;^'sgs for the generalised conductive flux 
in fluid of heat capacity cp and thermal conductivity x can be 
introduced: 

5"'™"''*=pcp(;r+A^sgs)Vr. (81) 

For the generalised convective flux y(<:°"'i)^ a closure migh t 
be based upon the super-adiabatic gradient (ICaniitnllTggi) . 
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fig.diff .png fig_diff_closr_vec .png 

(a) Explicit (b) = 0.065 (vectorial) 



f ig.dif f_locl_nb . png £ig_di£f _closr_scl . png 

(c) Q = C^' * (d) C, = 0.422 (scalar) 

Fig. 4. Turbulence energy isosurfaces as in Fig.[2with contour sections of the dimensionless flux magnitude of turbulent transport. 



Moreover, in some combustion problems or in simulations of 
multi-phase media the turbulent mixing of particle species is 
yet another challenge. These problems are left for future work. 

5. Turbulent burning in a box 

As a simple testing scenario, we performed LES of turbulent 
thermonuclear deflagration in degenerate carbon and oxygen. 



In these simulations, we utilised a greatly simplified reaction 
scheme, where the products of thermonuclear fusion are nickel 
and alpha particles. The thermonuclear burning zones prop- 
agate in a fashion similar to premixed chemical flames. For 
the chosen mass density, po « 2. 9 • lO^gcm'^, the width o f 
the flames is 6f ~ 0.006 cm (cf. iTimmes & WooslevliToQl . 
Hence, the flame fronts are appropriately represented by dis- 
continuities for the numerical resolution A = 2 ■ 10-' cm in the 
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simulations we run. The front propagation is numerically im- 
plement ed by mea ns of the level set method (Osher & Sethian 
Il988: Reinecke et al. 1999) . The domain of the flow is cubic 
with periodic boundary conditions (BCs). In this scenario, the 
burning process consumes all nuclear fuel within finite time. 
We set X = 216A = 4.32 km for the size of the domain, which 
is comparable to the resolution of the large scale supernova 
simulations to be discussed in paper 11. Since self-gravity is 
insignificant on length scales of the order of a few kilometres, 
we apply an external solenoidal force field in order to produce 
turbulent flow. Each Fourier mode of the force field is evolved 
as a distinct stochastic process of the Ornstein-Uhlenbeck type. 
The characteristic wavelength L of the forcing modes is half the 
size of the domain. L can be interpreted as integral length scale 
of the flow. An detailed description of the methodology and a 
discussion of numerous simulations is given in Schmid t et alJ 
ll2005ah . 

The LES of turbulent combustion is a particularly appro- 
priate case study for the performance of subgrid scale mod- 
els because the evolution of the system is strongly coupled to 
the SGS turbulence energy via the turbulent flame speed re- 
lation. For the notion of a turbulent fla me s peed seelPocheaul 
(ll9?4'), 'Nie mever & Hi llebrandt (199?) and Peters' ('l999'). In 
the framework of the filtering formalism, the underlying hy- 
pothesis is the following: If the flow is smoothed on a certain 
length scale A, then the effective propagation speed iturb(^) of 
a burning front is of the order of the turbulent velocity fluc- 
tuations v' ~ k^J^^^, provided that A » Zg- The length scale 
Ic is called the Gibson scale. It specifies the minimal size of 
turbulent eddies affecting the flame front propagation. In the 
context of a LES, we have ituA ~ ?sgs for the turbulent flame 
speed. Consequently, the SGS model determines the propaga- 
tion speed of turbulent flames. If Ic < A, on the other hand, 
the front propagation is determined by the microscopic con- 
ductivity of the fuel. The corresponding propagation speed is 
called the laminar flame speed and is denoted by si^^. Since 
iiam is determined by the balance between thermal conduction 
and thermonuclear heat generation, conduction effects are im- 
plicitly treated by the level set method. For this reason, we do 
not include the conduction terms in equation (I42> for the total 
energy. 

Both limiting cases of turbulent and laminar burning, re- 
spectivel y, are accommo dated in the flame speed relation pro- 
posed by[ 



i ', are accqmmo a 
Pocheaullll994l) : 



•Sturb 



1/2 



(82) 



The coefficient Ct is of the order unity and deter mines the ratio 
of Stuib and ^sgs in the turbulent burning regime. 'Peters' ('1999') 
propos es Q = 4/3, while Q ~ 20/3 is suggested by Kim et al. 
1I1999I) . Here we set Ct - 1 corre sponding to the asymp- 
totic relation iturb — ?sgs assumed by|N lemever & Hillebrandt" 
( ^1995). For a study of the influence of Ct see Schmidt et al. 

The laminar flame speed for an initial mass density of 
2.90 ■ 10^ gem"-' is iiam ~ 9.78 ■ 10^ cms"'. Choosing a charac- 
teristic velocity V = L/T = lOO^iam, where T is the autocorre- 
lation time of the stochastic force driving the flow, the Gibson 



scale becomes la ~ 10"^L ~ 0.1cm. Note that the Gibson 
length is still large compared to the flame thickness. Therefore, 
the internal structure of the burning zones is not disturbed by 
turbulent velocity fluctu ations, i.e. the flamelet regime of turbu- 
lent combustion applies jPetersI 19991) . 

Running a LES with the parameters outlined above and set- 
ting eight small ignition spots on a numerical grid of = 216^ 
cells, the expectation was that the burning process would ini- 
tially proceed slowly, but as turbulence was developing due to 
the action of the driving force, <7sgs would eventually exceed 
the laminar flame speed and substantially accelerate the flame 
propagation. Indeed, this is what can be seen in Fig. |6l which 
shows plots of statistical quantities as functions of time. The 
corresponding flame evolution is illustrated in the sequence 
of three-dimensional visualisations in Fig.0and|8j where the 
colour shading indicates the contour sections of ^sgs in loga- 
rithmic scaling. Initially, the spherical blobs of burning mate- 
rial are expanding slowly and become gradually elongated and 
folded by the onsetting flow which is produced by the driving 
force. As the SGS turbulence velocity g'sgs exceeds the laminar 
burning speed iiam in an increasing volume of space, the spa- 
tially averaged rate of nuclear energy release, (Pnuc), is increas- 
ing rapidly (Fig. |6}. Eventually, (Pnuc) assumes a peak value 
at dimensionless time f = t/T » 1.8 which coincides with 
the maximum of turbulence energy. Subsequently, the flow ap- 
proaches statistical equilibrium between mechanical produc- 
tion and dissipation of kinetic energy. Thus, the greater part of 
the fuel is burned within one large-eddy turn-over time of the 
turbulent flow. This observation in combination with the tight 
correlation between the growth of the mean rate of nuclear en- 
ergy release and the SGS turbulence velocity verifies that the 
burning process is dominated by turbulence. 

As a further indicator for the reliability of the SGS model, 
we varied the resolution in a sequence of LES, while main- 
taining the physical parameters unaltered. The resulting global 
statistics is shown in Fig.|9l In particular, the time evolution of 
(^'nuc) appears to be quite robust with respect to the numerical 
resolution. The deviations which can be discerned in the height, 
width and location of the peak are mostly a consequence of the 
different flow realisations due to the random nature of the driv- 
ing force. Actually, even if we had used identical sequences of 
random numbers to compute the stochastic force field in each 
simulation, the dependence of the time steps on the numerical 
resolution nevertheless would have produced different discrete 
realisations. Thus, we initialised the random number sequences 
differently and restricted the resolution study to statistical com- 
parisons. The evolution of the mass-weighted SGS turbulence 
velocity which is plotted in panel (b) of Fig. |51 reveals that tur- 
bulence is developing slightly faster in the case - 192^. 
This can be attributed to a somewhat larger root mean square 
force field during the first large-eddy turn-over in this simu- 
lation. Consequently, the burning process proceeds systemati- 
cally faster. Note, however, that the level of SGS turbulence be- 
comes monotonically lower with increasing resolution for the 
almost stationary flow at time t - 3.0. The deviations for the 
LES with the lowest resolution (A^ = 120-'), on the other hand, 
are likely to be spurious. For this reason, it would appear that 
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(a) Thermonuclear burning 



(b) SGS turbulence 



Fig. 6. Evolution of statistical quantities in a LES of thermonuclear deflagration in a cubic domain subject to periodic boundary 
conditions with = 216^ numerical cells. In panel (a) the spatially averaged rate of nuclear energy generation in combination 
with the mass fractions of unprocessed material (carbon and oxygen), alpha particles and nickel are plotted. The mean as well as 
the standard deviation of the SGS turbulence velocity g'sgs are shown in panel (b). 



the minimal resolution for sufficient convergence has to be set 
in between A? = 120^ and N ^ 160\ 

This conclusion is also supported by the turbulence energy 
spectra plotted in Fig. [TO] We computed the normalised energy 
spectrum functions for the transversal modes of the velocity 
fields after two integral time scales have elapsed. Details of 
the computation of discrete spectrum functions are discussed 
in (Schmidt et al.l (l2005bl) . One can clearly discern maxima 
in the vicinity of the normalised characteristic wave number 
k - Lkjln - 1.0 of the driving force. For the LES with 
greater than 120^, an inertial subrange emerges in the inter- 
val 2 < ^ < 6. The dimensionless cutoff wave number in the 
case — 216^ is ^ = 54. As demonstrated in Schmidt et al. 
ll2005bl) . the numerical dissipation of PPM, which was used to 
solve the hydrodynamical equations, noticeably smoothes the 
flow for wavenumbers k > 54/9 = 6. This is exactly what is 
observed in Fig.^| For A^ = 120^, on the other hand, virtu- 
ally all wavenumbers not directly affected by stochastic forc- 
ing are subject to numerical dissipation, i.e. there is no inertial 
subrange at all. Considering the more common power-of-two 
numbers of cells, a grid of N - 128^ cells will provide only 
marginally sufficient resolution, whereas one will be on the safe 
side with A^ = 256^. In paper II, however, it is shown that still 
higher resolutions might be required for LES of non-stationary 
inhomogeneous turbulence such as in the case of thermonu- 
clear supernova simulations. 

It is also argued in ISchmidt et all ll2005b') that the intrin- 
sic mean rate of dissipation produced by PPM closely agrees 
with the prediction of the Smagorinsky model for stationary 



isotropic turbulence. This suggests that the numerical dissipa- 
tion can be utilised as an implicit SGS model with regard to 
the velocity field. In fact, the LES presented in Fig.|6lEland|8] 
was computed without including the SGS stress term in the 
dynamical equation il6\ . while the total energy etot, which is 
conserved by PPM, was coupled to the SGS turbulence energy 
^sgs- One can think of k^^s as a buffer between the resolved ki- 
netic energy ^I^P and the internal energy eim. Apart from the 
energy budget, the SGS model influences the resolved dynam- 
ics via the turbulent flame speed. For the LES with varying 
resolution (Fig.lQlandllOL on the other hand, we applied com- 
plete coupling of the SGS model, i.e. the turbulent stress term 
in the momentum equation was included as well. Comparing 
Fig.|6l(a) and|5](a) for A^ - 216^, it appears that the burning 
process is slightly delayed in the latter case. As is discussed 
at length in Schmid t et alJ (l2005ah . the discrepancy can be at- 
tributed to a difficulty related to inverse energy transfer. Since 
backscattering injects energy on the smallest resolved scales, 
which are sizeably affected by numerical dissipation, the ki- 
netic energy added to the flow is more or less instantaneously 
converted into internal energy. Thus, the backscattering of en- 
ergy from subgrid scales to the resolved flow results in an arti- 
ficially enhanced dissipation which depletes turbulence energy. 
Using partial coupling, this unwanted effect is simply ignored. 
For consistency, one must then introduce a cutoff^ for the eddy- 
viscosity parameter Cy in order to dispose of negative viscosi- 
ties. Mending the shortcoming of the treatment of inverse en- 
ergy transfer is the subject of ongoing research. For the time be- 
ing, the partial coupling of the SGS model with backscattering 
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(a) t = 0.25T (b) t = 0.50T 
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Fig. 7. LES of thermonuclear deflagration in a cubic domain with periodic BCs. Shown are snapshots of the flame fronts with 
contour sections of the SGS turbulence velocity in logarithmic scahng. 



suppressed serves as a pragmatic solution in hydrodynamical 
simulations with PPM. 

6. Conclusion 

The localised SGS turbulence energy model olfers robustness 
and flexibiUty at relatively low computational cost. For this 



reason, it is particularly suitable for the application in LES 
of astrophysical fluid dynamics. The energy transfer from re- 
solved toward subgrid scales is modelled with the standard 
eddy-viscosity closure, where the closure parameter is com- 
puted from local properties of the flow. Hence, there are no a 
priori assumption about the resolved flow incorporated in the 
model. Non-local transport is treated with the down-gradient 
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Fig. 8. Fig.0continued. 



closure, using a constant statistical parameter. With a turbu- 
lent kinetic Prandtl number significantly larger than unity, it 
is possible to reproduce the magnitude of diffusive flux quite 
well. The rate of viscous dissipation appears to be particularly 
challenging. We found that a semi-statistical approach yields 
satisfactory results. 

The SGS model was implemented in a code for the LES 
of turbulent thermonuclear combustion in a periodic box using 



the piece-wise parabolic method (PPM) for the resolved hydro- 
dynamics and the level set method for the flame front propaga- 
tion. Since PPM produces significant numerical dissipation, we 
found it favourable to decouple the SGS model form the mo- 
mentum equation and suppressing inverse energy transfer from 
unresolved toward resolved scales. In this kind of application, 
the SGS turbulence energy serves as a buffer between the re- 
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(a) Rate of burning (b) SGS turbulence 

Fig. 9. Evolution of the mean dimensionless rate of nuclear energy generation (a) and the ratio of the mass-weighted mean SGS 
turbulence velocity to laminar burning speed (b) in a sequence of LES with varying resolution. 



solved kinetic energy and the internal energy and supplies a 
velocity scale for calculation of the turbulent burning speed. 

Furthermore, gravitational and thermal effects can be in- 
cluded in the SGS model, although closures specific to a cer- 
tain physical system have to be formulated. An example is pre- 
sented in paper II, where the application of the SGS model to 
Rayleigh-Taylor-driven thermonuclear combustion in type la 
supernova is discussed. Adapting the model to other applica- 
tions, possibly with different numerical techniques, is the goal 
of on-going research. 
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